/**

This file is part of MaCI/GIMnet.

MaCI/GIMnet is free software: you can redistribute it and/or modify it 
under the terms of the GNU Lesser General Public License as published 
by the Free Software Foundation, either version 3 of the License, or 
(at your option) any later version.

MaCI/GIMnet is distributed in the hope that it will be useful, but WITHOUT 
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
License for more details.

You should have received a copy of the GNU Lesser General Public 
License along with GIMnet. (See COPYING.LESSER) If not, see 
<http://www.gnu.org/licenses/>.

**/
#include <stdlib.h>
#include <math.h>
#include <mex.h>


#define zeros(X,Y) mxCreateDoubleMatrix(X,Y,mxREAL);
#define sendup(X) mxSetName(X,#X); mexPutArray(X,"base");
#define MAXF 1000
#define CONV(a,b) (abs(a-b))
#define BASETYPE2 double
#define BASETYPE  unsigned char
#define PI 3.14159265358979

/*mex -v kmean.c */
void mgem1(double *x,double *w, double *p,double *lam, double *W, double *m, double *mu, int N, int numpoints, int C)
{
    double e,sum,cumsum;
    int i,j,k;
    
    //calculate likelihood of each point
    //w is numcenters x numpoints matrix
    //P(c|x,mu)=P(x|c,mu)P(c)/P(x|c)
    for(i=0;i<numpoints;i++)
    {
        lam[i]=0;
        for(j=0;j<C;j++)
        {
            e=x[i*N]-m[j*N];
            sum=e*e;
            for(k=1;k<N;k++)
            {
                e=x[i*N+k]-m[j*N+k];
                sum+=e*e;
            }
            sum=sum*p[j];
            lam[i]+=sum;
            w[i*C+j]=sum;
        }
        cumsum+=lam[i];
        //normalize likelihoods across rows
        for(j=0;j<C;j++)
        {
            w[i*C+j]=exp(w[i*C+j]-lam[i]);
        }
    }
    //normalize likelihoods across columns
    for(j=0;j<C;j++)
    {
        W[j]=w[j];
        for(i=1;i<numpoints;i++)
        {
            W[j]+=w[i*C+j];
        }
    }
    //recalculate means
    for(i=0;i<numpoints;i++)
    {
        for(j=0;j<C;j++)
        {
            for(k=0;k<N;k++)
            {
                mu[j*N+k]+=w[i*C+j]*x[i*N+k];
                
            }
        }
    }
    
}

void kmeans(double *x, double *mu, double *r,double *R,double *res, int N, int K, int numPoints, double *iter)
{
    double e,sum,min,fac;
    int i,j,k,a,conv=0,cyc=0;

    while(conv==0)
    {
        conv=1;
        //assigment step
        for(k=0;k<K;k++)
        {
            R[k]=0;
            res[k]=0;
        }
        for(i=0;i<numPoints;i++)
        {
            e=x[i*N]-mu[0];
            min=e*e;
            for(j=1;j<N;j++)
            {
                e=x[i*N+j]-mu[j];
                min+=(e*e);
            }
            a=0;
            for(k=1;k<K;k++)
            {
                e=x[i*N]-mu[k*N];
                sum=(e*e);
                for(j=1;j<N;j++)
                {
                    e=x[i*N+j]-mu[k*N+j];
                    sum+=(e*e);
                }
       
                if(sum<min)
                {
                    min=sum;
                    a=k;
                }
            }

            
            R[a]++;
            res[a]+=min;
            if(r[i]!=a)
                conv=0;
            r[i]=a;
        }
            //update
            
        
        for(k=0;k<K;k++)
        {
            for(j=0;j<N;j++)
            {
                mu[k*N+j]=0;
            }
        }

        for(i=0;i<numPoints;i++)
        {
            a=(int)r[i];
            for(j=0;j<N;j++)
            {
                mu[a*N+j]+=x[i*N+j];
            }
            
        }
        for(k=0;k<K;k++)
        {
            fac=1./R[k];
            for(j=0;j<N;j++)
            {
                mu[k*N+j]=mu[k*N+j]*fac;
            }
        }
        cyc++;
        (*iter)=(double)cyc;
        if(cyc>40)
            break;
       // conv=1;
    }
}
/******************************************************/
void mexFunction(
int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]
)
{  

  
  const mxArray     *mmu, *mp, *mx;
  mxArray  *ml,*mmuh,*mcyc,*mres,*mR;
 
  double *x,*p,*l, *mu,*muh,*tmp,*cyc,*res,*R;
 
  int    i, N, K, numPoints;
  
  

  
  static int	expectRhs = 3;
  static int	expectLhs = 5;
  char		buffer[100];

  /* Check for proper number of arguments */

  if (nrhs != expectRhs) 
  {
    mexPrintf("kmeans requires %d input arguments. Got %d",
	    expectRhs, nrhs);
  } else if (nlhs > expectLhs) 
  {
    mexPrintf("kmeans requires %d outputs arguments. Got %d",
	    expectLhs, nlhs);
  }
  
  i=0;
  mmu = prhs[i++];
  mp = prhs[i++];
  mx = prhs[i++];
 // cyc = mxGetScalar(prhs[i++]);




 
  K = mxGetN(mmu); /* colums */
  N = mxGetM(mmu); /* rows */


  numPoints = mxGetN(mx); 
  mexPrintf("num:%d\n",numPoints);
  
  mcyc=zeros(2,1);
  cyc=mxGetPr(mcyc);
  ml=zeros(numPoints,1);
  l=mxGetPr(ml);
  mmuh=zeros(N,K);
  mres=zeros(K,1);
  mR=zeros(K,1);
  res=mxGetPr(mres);
  R=mxGetPr(mR);
  muh=mxGetPr(mmuh);
  mu=mxGetPr(mmu); 
  p=mxGetPr(mp); 
  x=mxGetPr(mx); 

  
 
  for(i=0;i<N*K;i++)
      muh[i]=mu[i];
  for(i=0;i<numPoint2s;i++)
      l[i]=p[i];
/**
*  K = Number of klusters 
*  N = Number of dimension 
*  NumPoints = Number of data points

* @p  *x  The data (N * Numpoints)[x x x x x; y y y y ; a a a a]
* @p  *mu Average points  (N * K) 
* @p  *r  (NumPoints) = Label mihin klusteriin kuuluu (Klusterin indeksi jokaiselle datapisteelle)
* @p  *R  K*1 Kuinka monta pointsia on klustereissa
* @p  *res residuaali = neliösumma = pisteiden etäisyys omasta klusterista (K)
* @p  N
* @p K
* @p numPoints
* @p *iter Number of iterations
* 
* 
**/
	void kmeans(double *x, double *mu, double *r,double *R,double *res, int N, int K, int numPoints, double *iter)
  kmeans(x,muh, l,R,res, N, K, numPoints, cyc);
/*    for(i=0;i<cyc;i++)
    {
        sum=0;
        for(k=0;k<M;k++)
        {
            lC[k]=C[N*k];
            for(l=1;l<N;l++)
                lC[k]*=C[N*k+l];
        }
        for(j=0;j<numPoints;j++)
        {
        // calculate weights wi=exp(-1/2*(pts-mui)'*iCi*(pts-mui))
            
            e=pts[j*N]*pts[j*N]*iC[0];
            for(k=1;k>N;k++)
                e+=pts[j*N+k]*pts[j*N+k]*iC[k];
            sum+=e;
        // calculate new means from weighted datapoints
        // calculate new covariances from weighted datapoints
        }
    }
*/
  

  i = 0;
  

  plhs[i++] = mmuh;
  plhs[i++] = ml;
  plhs[i++] = mcyc;
  plhs[i++] = mres;
  plhs[i++] = mR;



 
}
